package scu.maqiang.electromagnetics.waveguide;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class RectangleWaveguide {
    public static void main(String[] args) {
        //RectangleWaveguide();
        coaxialLine();
    }

    public static void RectangleWaveguide() {
        Mesh2T3 mesh = new Mesh2T3().square2D(20, 10, T3Type.Left);
        mesh.scale(2.0, 1.0);

        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SRMatrix M = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, M);
        double[] x = MVO.rand(fs.GetNdof());

        int nVec = 20;
        double[] D = new double[nVec];
        double[][] V = new double[nVec][fs.GetNdof()];
        double[][] dxyV = new double[2][fs.GetNdof()];
        fs.applyBC_MBN(A, x, Direct.All, 0.0, 1, 2, 3, 4);
        EigenSSolver solver = new EigenSSolver();
        solver.inversePowerTran_DSS(A, M, x, 10.0, nVec, D, V);
        Arrays.setAll(D, i -> 2 * Math.sqrt(D[i]));
        System.out.println(MVO.toString(D));
        double[][] result = new double[3][fs.GetNdof()];
        result[0] = V[0];
        fs.computeGradient(result[0], dxyV);
        Arrays.setAll(result[1], i-> -dxyV[0][i]);
        Arrays.setAll(result[2], i-> -dxyV[1][i]);
        mesh.toTecplot("ElectricField.dat", result);
    }

    public static void coaxialLine() {
        Mesh2T3 mesh = new Mesh2T3().ring(1.0, 4.0, 76, 20);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        SRMatrix M = new SRMatrix(fs.GetNdof());
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.assembleMass(new double[]{1.0}, BVPType.CONSIST_MASS, M);
        double[] x = MVO.rand(fs.GetNdof());

        int nVec = 20;
        double[] D = new double[nVec];
        double[][] V = new double[nVec][fs.GetNdof()];
        double[][] dxyV = new double[2][fs.GetNdof()];
//        fs.applyBC_MBN(A, x, Direct.All, 0.0, 1, 2, 3, 4);
        EigenSSolver solver = new EigenSSolver();
        solver.inversePowerTran_DSS(A, M, x, 10.0, nVec, D, V);
        Arrays.setAll(D, i -> 1 * Math.sqrt(D[i]));
        System.out.println(MVO.toString(D));
        double[][] result = new double[3][fs.GetNdof()];
        result[0] = V[0];
        fs.computeGradient(result[0], dxyV);
        Arrays.setAll(result[1], i-> -dxyV[0][i]);
        Arrays.setAll(result[2], i-> -dxyV[1][i]);
        mesh.toTecplot("CoAxialLine.dat", result);
    }
}
